#additional stuff that referees asked for: expand "extra seats" in 2010.
function runExtraCounterfactuals(theta,dataset,counterfactuals,lv=dataset[2012][2],temps=dataset[2012][3])
    #set up counterfactual stuff here
    (cf_constants, cf_temps) = counterfactuals[2012]
    cf_temps.alphadraws[1:end-1,:] .= rand.()
    cc = dataset[2012][1]
    cap0 = [cf_constants.seats+cf_constants.sobrecupo; 0] #baseline, initial match
    II = size(lv.U,2)
    cap0_after = [cf_constants.seats; II]
    cap_efficient = [max.(cf_constants.seats, [sum(cc.enroll.==jj) for jj=1:size(lv.U,1)]); II]
    cap0_after_noG8 = copy(cap0_after)
    cap0_after_noG8[findall(cf_constants.G8)] .= cap_efficient[findall(cf_constants.G8)]
    df_individuals = DataFrame(ind=1:II, type=cc.typ)
    df_means = DataFrame(inds=["All","Male Private","Male Public","Female Private","Female Public"])
    df_byIter = Dict{Symbol,Array{Float64,1}}()
    dfs = (df_byIter,df_means,df_individuals)
    #
    platform = .!(cf_constants.G8)
    asif2011 = twoStageMatch(dfs,:asif2011,lv,cap0.*[platform;false],cap0_after_noG8,theta.alpha,platform,cf_constants,cf_temps,cc,theta)
    for extraSeats in 0.05 : 0.05 : 0.30
        eee = Int(round(100 * extraSeats))
        cap_pre = Int.(round.(cap0.*[platform;false].*(1 .+ extraSeats)))
        asif2011_e = twoStageMatch(dfs,Symbol("asif2011_addseats_$eee"),lv,cap_pre,cap0_after_noG8,theta.alpha,platform,cf_constants,cf_temps,cc,theta)
    end
    df_means, df_individuals, df_byIter
end


function runCounterfactuals(theta,dataset,counterfactuals,lv=dataset[2012][2],temps=dataset[2012][3]; dropInstitutions=false, runMinimalCounterfactuals=runMinimalCounterfactuals)
    #set up counterfactual stuff here
    (cf_constants, cf_temps) = counterfactuals[2012]
    cf_temps.alphadraws[1:end-1,:] .= rand.()
    cc = dataset[2012][1]
    cap0 = [cf_constants.seats+cf_constants.sobrecupo; 0] #baseline, initial match
    II = size(lv.U,2)
    cap0_after = [cf_constants.seats; II]
    cap_efficient = [max.(cf_constants.seats, [sum(cc.enroll.==jj) for jj=1:size(lv.U,1)]); II]
    cap0_after_noG8 = copy(cap0_after)
    cap0_after_noG8[findall(cf_constants.G8)] .= cap_efficient[findall(cf_constants.G8)]
    df_individuals = DataFrame(ind=1:II, type=cc.typ)
    df_means = DataFrame(inds=["All","Male Private","Male Public","Female Private","Female Public"])
    df_byIter = Dict{Symbol,Array{Float64,1}}()
    dfs = (df_byIter,df_means,df_individuals)
    #run them!
    #maximize utility (mess with priorities to do so)
    cfc = deepcopy(cf_constants)
    for jj=1:size(lv.U,1)
        uj_scaled = lv.U[jj,:] ./ abs.(theta.beta_ij[1,cc.typ])
        cfc.priorityEligible[:,jj] .= invperm(sortperm(uj_scaled)) .* (cf_constants.priorityEligible[:,jj] .> 0)
    end
    maxutility = oneStageMatch(dfs,:maxutility,lv,[cap_efficient;II],cfc,cf_temps,cc,theta)
    #third: everything else (simpler), condition on utilities in data.
    zeroalpha = twoStageMatch(dfs,:zeroalpha,lv,cap0,cap0_after,theta.alpha,cc.platform,cf_constants,cf_temps,cc,theta; scale_frictions=0.0)
    cap_efficient = [max(cf_constants.seats[jj], sum(zeroalpha.placement .==jj)) for jj=1:size(lv.U,1)]
    efficient = oneStageMatch(dfs,:efficient,lv,[cap_efficient;II],cf_constants,cf_temps,cc,theta)
    asInData = getAsInData(cc,lv,theta)
    appendResults!(dfs, asInData, Symbol("asInData"), cc, cf_constants)
    benchmark = twoStageMatch(dfs,:benchmark,lv,cap0,cap0_after,theta.alpha,cc.platform,cf_constants,cf_temps,cc,theta)
    platform = .!(cf_constants.G8)
    asif2011 = twoStageMatch(dfs,:asif2011,lv,cap0.*[platform;false],cap0_after_noG8,theta.alpha,platform,cf_constants,cf_temps,cc,theta)
    equalalpha = let
        alpha = copy(theta.alpha)
        alpha .= theta.alpha[:,1]
        twoStageMatch(dfs,:equalalpha,lv,cap0.*[platform;false],cap0_after_noG8,alpha,platform,cf_constants,cf_temps,cc,theta)
    end
    smallalpha = twoStageMatch(dfs,:smallalpha,lv,cap0.*[platform;false],cap0_after_noG8,theta.alpha,platform,cf_constants,cf_temps,cc,theta,scale_frictions=0.1)
    smallalpha2012 = twoStageMatch(dfs,:smallalpha2012,lv,cap0 ,cap0_after,theta.alpha,cc.platform,cf_constants,cf_temps,cc,theta,scale_frictions=0.1)
    #"free college for low SES"
    lv_free, lv_free_G25 = make_lv_freecollege(theta,cc,lv,cf_constants)
    freecollege = twoStageMatch(dfs,:freecollege_lowSES,lv_free,cap0,cap0_after,theta.alpha,cc.platform,cf_constants,cf_temps,cc,theta)
    freecollege2011 = twoStageMatch(dfs,:freecollege_lowSES_asif2011,lv_free_G25,cap0.*[platform;false],cap0_after_noG8,theta.alpha,platform,cf_constants,cf_temps,cc,theta)
    #sequence of shrinking alphas
    if !runMinimalCounterfactuals
        for scale = 0.0 : .1 : .9
            s = twoStageMatch(dfs,Symbol("scaleAlpha_$(Int(10*scale))"),lv,cap0 ,cap0_after,theta.alpha,cc.platform,cf_constants,cf_temps,cc,theta; scale_frictions=scale)
            s2011 = twoStageMatch(dfs,Symbol("scaleAlpha2011_$(Int(10*scale))"),lv,cap0.*[platform;false],cap0_after_noG8,theta.alpha,platform,cf_constants,cf_temps,cc,theta; scale_frictions=scale)
        end
        #new "drop programs by selectivity" cfctl
        for dec in 1:maximum(programselectivitybin)
            pl = (programselectivitybin .!= dec)
            c = twoStageMatch(dfs,Symbol("dropSelectivityBin$dec"),lv,cap0.*[pl;false],cap0_after,theta.alpha,pl,cf_constants,cf_temps,cc,theta; addIndividualOutcomes=false)
        end
    end
    #old "drop institutions" cfctl
    if dropInstitutions
        for jj in unique(cf_constants.institutionID)  #0 is "don't drop anything" cfctl
            println("counterfactual: excluding institutionID $jj")
            J = size(lv.U,1)
            keep_jj_off = [trues(J);false]
            keep_jj_off[findall(cf_constants.institutionID .== jj)] .= false
            result = twoStageMatch(dfs,Symbol("dropInstitution$jj"),lv,cap0.*keep_jj_off,cap0_after,theta.alpha,keep_jj_off[1:J],cf_constants,cf_temps,cc,theta; addIndividualOutcomes=false)
        end
    end
    df_means, df_individuals, df_byIter
end

function setupInitialDA!(cf_temps,cf_constants,lv,platform,cc,theta)
    (J,I) = size(lv.U)
    #fill in Ufull, students' prefs
    Threads.@threads for ii=1:I
        typei = cc.typ[ii]
        Ui = view(cf_temps.Ufull,:,ii) #assume comes from 2012 when everything is on-platform
        @simd for jj=1:J
            Ui[jj] = lv.U[jj,ii]
            # if !platform[jj]
            #     Ui[jj] += theta.beta_fixed[1,typei] #adjust utility of off-platform programs
            # end
        end
        Ui[J+1] = lv.U0[2,ii]
    end
    #fill in initial best offer and utility from it w/ first oo
    cf_temps.point_at .= 0
    cf_temps.U_pointed .= view(lv.U0,1,:)
    cf_temps.enroll .= 0
    cf_temps.programPrefs .= cf_constants.priorityEligible
    cf_temps.proposalcount .= 0
    getProgramROL!(cf_temps)
end

function dropPeopleWhoDidNotApply!(cf_temps)
    (II,Jtilde) = size(cf_temps.programROL)
    @assert 1000000 > maximum(cf_temps.programPrefs)
    cutoff = 1000000 .* ones(eltype(cf_temps.ProgramPrefs),Jtilde)
    for jj=1:Jtilde
        if cf_temps.proposalcount[jj] > 0
            lastAdmitted = cf_temps.programROL[cf_temps.proposalcount[jj],jj]
            cutoff[jj] = cf_temps.programPrefs[lastAdmitted,jj]
        end
    end
    cf_temps.programPrefs[cf_temps.Ufull' .< cf_temps.U_pointed] .= 0
    getProgramROL!(cf_temps)
    for jj=1:Jtilde
        if cf_temps.proposalcount[jj] > 0
            cf_temps.proposalcount[jj] = 0
            for npropose = 1:II
                ind_propose = cf_temps.programROL[npropose,jj]
                if (ind_propose==0) || (cf_temps.programPrefs[ind_propose,jj] < cutoff[jj])
                    cf_temps.proposalcount[jj] = npropose-1
                    break
                end
                ind_propose==II && (cf_temps.proposalcount[jj] = II)
            end
        end
    end
end

function getProgramROL!(cf_temps)
    (Jtilde,I) = size(cf_temps.Ufull)
    Threads.@threads for jj=1:Jtilde
        _getApp!(view(cf_temps.programROL,:,jj), view(cf_temps.programPrefs,:,jj), 0)
    end
end

function _getApp!(app,Ui,u00)
    Jtilde = length(Ui)
    rank_oo0 = 1
    @simd for t=1:Jtilde
        rank_oo0 += (Ui[t] > u00)
    end
    partialsortperm!(app,Ui,1:rank_oo0-1,rev=true)
    @simd for t = rank_oo0:Jtilde
        app[t] = 0
    end
end

#student-proposing DA
function runSPDA!(cf_constants,cf_temps,cap,maxiter=500000)
    II, Jtilde = size(cf_constants.priorityEligible)
    @assert Jtilde == length(cap)
    # construct student ROLs.  [rr,ii] cell is ii's rrth preference.
    Uoo = cf_temps.U_pointed  #initially contains first outside option!
    studentROL = zeros(Int,Jtilde,II)
    Threads.@threads for ii=1:II
        sortperm!(view(studentROL,:,ii),view(cf_temps.Ufull,:,ii),rev=true)
    end
    # run it!
    bar = ones(Jtilde) #minimum score needed for admission
    bar_new = ones(Jtilde)
    for iter=1:maxiter
        println("SPDA iter $iter")
        #first, students apply
        cf_temps.proposalcount .= 0
        cf_temps.point_at .= 0
        Threads.@threads for ii=1:II
            for rr=1:Jtilde
                jj = studentROL[rr,ii]
                cf_temps.Ufull[jj,ii] <= Uoo[ii] && break
                if cf_constants.priorityEligible[ii,jj] >= bar[jj]
                    cf_temps.point_at[ii] = jj
                    break
                end
            end
        end
        #second, programs reject
        Threads.@threads for jj=1:Jtilde
            applicants = findall(cf_temps.point_at .== jj)
            pri = cf_constants.priorityEligible[applicants,jj]
            if length(pri) >= cap[jj]
                if cap[jj] == 0
                    bar_new[jj] = Inf
                else
                    bar_new[jj] = partialsort(pri,cap[jj],rev=true)
                end
            end
        end
        @assert all(bar_new .>= bar)
        converged = all(bar_new .== bar)
        bar .= bar_new
        converged && return true
    end
    return false
end

#college-proposing DA
function runDA!(cf_temps,capacity,maxiter=500000)
    for iter=1:maxiter
        updateEnrollment!(cf_temps)
        converged = collegesPropose!(cf_temps,capacity)
        converged && return true
    end
    return false
end

function updateEnrollment!(cf_temps)
    cf_temps.enroll .= 0
    for jj in cf_temps.point_at
        jj>0 && (cf_temps.enroll[jj] += 1)
    end
end

function collegesPropose!(cf_temps,capacity)
    converged = true
    (Jtilde,II) = size(cf_temps.Ufull)
    proposalcount = cf_temps.proposalcount #rank of last offer made by j
    enrollment = cf_temps.enroll
    for jj=1:Jtilde
        programROL = cf_temps.programROL
        if enrollment[jj] < capacity[jj] && proposalcount[jj] < II && programROL[proposalcount[jj]+1,jj] != 0
            stop_at = proposalcount[jj] + capacity[jj] - enrollment[jj]
            if stop_at < II && programROL[stop_at,jj] > 0
                while programROL[stop_at+1,jj]==programROL[stop_at,jj]
                    stop_at += 1 #deal with ties: everyone who ties for last offer gets one
                end
            end
            for ind = proposalcount[jj]+1:stop_at
                ind > II && break
                ii = programROL[ind,jj]
                ii == 0 && break #exhausted all acceptable applicants
                converged = false
                proposalcount[jj] += 1
                Uij = cf_temps.Ufull[jj,ii]
                if Uij > cf_temps.U_pointed[ii]
                    cf_temps.U_pointed[ii] = Uij
                    cf_temps.point_at[ii] = jj
                end
            end
        end
    end
    converged
end

#implement aftermarket frictions by making matches with bad alpha draws unattractive to students
function frictions!(cc,cf_temps,alpha,platform,typelist,scaling=1.0)
    (Jtilde, I) = size(cf_temps.Ufull)
    J = length(platform)
    @assert J+1==Jtilde
    prNotAvailable = zeros(J)
    Xfr = zeros(J,4)
    for ii=1:I
        typ = typelist[ii]
        #Xfr .= cc.Xfriction_full[ii]
        for j=1:J
            if cc.G8[j]
                Xfr[j,1] = 0.0
                Xfr[j,2] = platform[j]
                Xfr[j,3] = !platform[j]
            else
                Xfr[j,1] = 1.0
                Xfr[j,2] = 0.0
                Xfr[j,3] = 0.0
            end
            Xfr[j,4] = cc.Xfriction_full[ii][j,4]
        end
        mul!(prNotAvailable, Xfr, view(alpha,:,typ)) #index
        prNotAvailable .= cdf.(Normal.(prNotAvailable,1),0) #pr( index + N(0,1) < 0)
        prNotAvailable .*= scaling
        for jj=1:J
            (cf_temps.point_at[ii] != jj) && (cf_temps.alphadraws[jj,ii] < prNotAvailable[jj]) && (cf_temps.Ufull[jj,ii] = cf_temps.U_pointed[ii] - 10000)
        end
    end
end

function getWelfare(cc,lv,cf_temps,theta)
    u0 = vec( maximum(lv.U0,dims=1) )
    welfare = (cf_temps.U_pointed .- u0) ./ abs.(theta.beta_ij[1,cc.typ])
end

#faster in-place version
function getWelfare!(u0,welfare,cc,lv,cf_temps,theta)
    u0 .= max.(view(lv.U0,1,:), view(lv.U0,2,:))
    welfare .= (cf_temps.U_pointed .- u0)
    for ii=1:length(welfare)
        welfare[ii] /= abs(theta.beta_ij[1,cc.typ[ii]])
    end
    welfare
end

function twoStageMatch(dfs,name,
        lv,cap_initial,cap_after,alpha,platform,
        cf_constants=cf_constants,cf_temps=cf_temps,cc=cc,theta=theta;
            maxIterations=60000, addIndividualOutcomes=true, scale_frictions=1.0, storeInitialPlacement=false)
    setupInitialDA!(cf_temps,cf_constants,lv,platform,cc,theta)
    dropPeopleWhoDidNotApply!(cf_temps)
    runDA!(cf_temps,cap_initial)
    placement0 = copy(cf_temps.point_at)
    n_renege = sum( (cf_temps.U_pointed .< cf_temps.Ufull[end,:]) .& (cf_temps.point_at .> 0))
    #dropPeopleWhoDidNotApply!(cf_temps)
    frictions!(cc,cf_temps,alpha,platform,cc.typ,scale_frictions)
    welfareByIter = Float64[]
    sharePlacedByIter = Float64[]
    shareImprovedByIter = Int[]
    u0 = zeros(size(lv.U,2))
    welf = zeros(size(lv.U,2))
    welf_old = (cf_temps.U_pointed .- vec(maximum(lv.U0,dims=1))) ./ abs.(theta.beta_ij[1,cc.typ])
    for t=1:maxIterations::Int
        converged = runDA!(cf_temps,cap_after,1)
        getWelfare!(u0,welf,cc,lv,cf_temps,theta)
        push!(welfareByIter,mean(welf))
        sp=0 #n placed
        si=0 #n improved
        J = size(lv.U,1)
        for (ii,p) in enumerate(cf_temps.point_at)
            sp += (0 < p < J+1)
            #si += (0 < p < J+1) && (welf[ii]>welf_old[ii])
        end
        push!(sharePlacedByIter,sp/length(welf))
        #push!(shareImprovedByIter,si)
        #welf_old .= welf
        converged && break
    end
    placement = copy(cf_temps.point_at)
    Jtilde = size(cf_temps.Ufull,1)
    placement[placement .== Jtilde] .= 0
    welfare = getWelfare(cc,lv,cf_temps,theta)
    graduation = simulateOutcome(cc,lv,placement,theta,platform)
    println("final DA: mean welfare = $(mean(welfare)), share new placement = $(mean(placement0 .!= placement)), grad = $(sum(graduation)/sum(placement .> 0))")
    cfctl = (placement0=placement0, placement=placement, welfare=welfare, graduation=graduation,
            welfareByIter=welfareByIter,sharePlacedByIter=sharePlacedByIter,
            shareImprovedByIter=shareImprovedByIter)
    appendResults!(dfs, cfctl, name, cc, cf_constants, addIndividualOutcomes, storeInitialPlacement)
    return cfctl
end

function oneStageMatch(dfs,name,lv,cap,cf_constants=cf_constants,cf_temps=cf_temps,cc=cc,theta=theta)
    Jtilde = size(cf_temps.Ufull,1)
    platform = trues(Jtilde-1)
    setupInitialDA!(cf_temps,cf_constants,lv,platform,cc,theta)
    dropPeopleWhoDidNotApply!(cf_temps)
    runDA!(cf_temps,cap)
    placement = copy(cf_temps.point_at)
    placement[placement .== Jtilde] .= 0
    welfare = getWelfare(cc,lv,cf_temps,theta)
    graduation = simulateOutcome(cc,lv,placement,theta,platform)
    println("One-stage DA: mean welfare = $(mean(welfare)).")
    cfctl = (placement=placement,welfare=welfare,graduation=graduation)
    appendResults!(dfs, cfctl, name, cc, cf_constants, true)
    return cfctl
end

function getAsInData(cc,lv,theta)
    (J,I) = size(lv.U)
    enroll = copy(cc.enroll)
    welfare = zeros(I)
    for ii=1:I
        if enroll[ii] == 0
            welfare[ii] = 0 #max(lv.U0[1,ii],lv.U0[2,ii])
        else
            u0 = max(lv.U0[1,ii],lv.U0[2,ii])
            welfare[ii] = (lv.U[enroll[ii],ii]-u0)/theta.beta_ij[1,cc.typ[ii]]
        end
    end
    enroll[enroll.==J+1] .= 0
    graduation = 1.0 .* cc.Grad6
    (placement=enroll,welfare=welfare,graduation=graduation)
end

function summarizeWelfare(cc,uvec,enrollvec,graduationvec)
    #betterenroll_typ = zeros(cc.ntypes)
    graduation_typ = zeros(cc.ntypes)
    enroll_typ = zeros(cc.ntypes)
    welf_typ = zeros(cc.ntypes)
    n_typ = zeros(Int,cc.ntypes)
    J = size(cc.J_ub,1)
    for (u,e,g,typ) in zip(uvec, enrollvec, graduationvec, cc.typ)
        enroll_typ[typ] += (0 < e < J+1)
        welf_typ[typ] += u
        graduation_typ[typ] += g
        n_typ[typ] += 1
    end
    welfare = [sum(welf_typ)./sum(n_typ); welf_typ ./ n_typ]
    enroll = [sum(enroll_typ)./sum(n_typ); enroll_typ ./ n_typ]
    grad = [sum(graduation_typ)./sum(enroll_typ); graduation_typ ./ enroll_typ]
    return welfare, enroll, grad
end

#appendResults!(dfs::Nothing,cfctl,name, cc, cf_constants, addIndividualOutcomes=true) = nothing
function appendResults!(dfs, cfctl, name, cc, cf_constants, addIndividualOutcomes=true, storeInitialPlacement=false)
    println("summarizing results from counterfactual: $name")
    (df_byIter, df_means, df_individuals) = dfs
    mean_welf, mean_enroll, mean_grad = summarizeWelfare(cc,cfctl.welfare,cfctl.placement,cfctl.graduation)
    df_means[!,Symbol("$(name)_welfare")] = mean_welf
    df_means[!,Symbol("$(name)_enroll")] = mean_enroll
    df_means[!,Symbol("$(name)_graduation")] = mean_grad
    if haskey(cfctl,:welfareByIter)
        df_byIter[Symbol("$(name)_welfare")] = cfctl.welfareByIter
        df_byIter[Symbol("$(name)_enroll")] = cfctl.sharePlacedByIter
    else
        df_byIter[Symbol("$(name)_welfare")] = [mean_welf[1],]
        df_byIter[Symbol("$(name)_enroll")] = [mean_enroll[1],]
    end
    if addIndividualOutcomes
        if storeInitialPlacement
            df_individuals[!,Symbol("$(name)_placement")] = cfctl.placement0
        end
        df_individuals[!,Symbol("$(name)_enroll")] = cfctl.placement
        df_individuals[!,Symbol("$(name)_welfare")] = cfctl.welfare
        df_individuals[!,Symbol("$(name)_graduate")] = cfctl.graduation
    end
end

function simulateOutcome(cc,lv,placement,theta,platform)
    (J,II) = size(lv.U)
	nmatch = size(cc.Xij[1],2)
	nfix = size(cc.Xfixed,2)
	noo = size(cc.Xoo,2)
    Xoutcome = zeros(1+nmatch+nfix+noo)
    grad6 = zeros(II)
    for ii=1:II
        e=placement[ii]
        typei = cc.typ[ii]
        if e > 0
            #Xoutcome[1] = log(lv.U0[1,ii])
            Xoutcome[1] = lv.U[e,ii]
            Xoutcome[2:noo+1] .= view(cc.Xoo,ii,:)
            Xoutcome[noo+2: noo+nmatch+1] .= view(cc.Xij[ii],e,:)
            Xoutcome[noo+nmatch+2:end] .= view(cc.Xfixed,e,:)
            #Xoutcome[noo+nmatch+3] = 1-platform[e] #1st col of Xfixed is 1(offplatform)
            beta = view(theta.betaOutcome,:,typei)
            xb = dot(Xoutcome,beta)
            grad6[ii] = ccdf(Normal(xb,1.0),0.0)
        end
    end
    grad6
end

function make_lv_freecollege(theta,cc,lv,cf_constants)
    lv_free = (U=copy(lv.U), U0=copy(lv.U0), Astar=copy(lv.Astar), expenditure=zeros(size(lv.U)))
    lv_free_G25 = (U=copy(lv.U), U0=copy(lv.U0), Astar=copy(lv.Astar), expenditure=zeros(size(lv.U)))
    for ii=1:size(lv_free.U,2)
        for jj=1:size(lv.U,1)
            netprice = cf_constants.Pij[ii,jj] #cc.Xij[ii][jj,1]
            u_gain = -netprice*theta.beta_ij[1,cc.typ[ii]] #+ cc.Xij[ii][jj,2]*theta.beta_ij[2,cc.typ[ii]]
            lv_free.U[jj,ii] += u_gain
            lv_free.expenditure[jj,ii] = netprice
            if !cf_constants.G8[jj]
                lv_free_G25.U[jj,ii] += u_gain
                lv_free_G25.expenditure[jj,ii] = lv_free.expenditure[jj,ii]
            end
        end
    end
    return lv_free, lv_free_G25
end
